ESC-TR-2003-076 


Project  Report 
PCA-IRT-3 

Preliminary  Design  Review:  GMTI 
Narrowband  for  the  Basic  PCA  Integrated 

Radar-Tracker  Application 


A. I.  Reuther 


27  February  2003 
Issued  6  February  2004 


Lincoln  Laboratory 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
Lexington,  Massachusetts 


Prepared  for  the  Defense  Advanced  Research  Projects  Agency 
under  Air  Force  Contract  F19628-00-C-0002. 


Approved  for  public  release;  distribution  is  unlimited. 


This  report  is  baser!  on  sUidies  performed  at  Lincoln  Laboratory,  a  center  for 
research  operated  by  Massachusetts  Institute  of  Technology.  This  -work  was 
sponsoi'ed  by  DARPA/ITO  under  Air  Force  Contract  F19628-00-C-0002. 

Tliis  report  may  be  reproduced  to  satisfy  needs  of  U.S.  Government  agencies. 


The  ESC  Public  Affairs  OflSce  has  reviewed  this  report,  and 
itis  releasable  to  the  National  Technical  Information  Service, 
where  it  wiU  be  available  to  the  general  public,  including 
foreign  nationals. 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 
FOR  THE  C0MRL4NDER 


"^0 

Cary  Adtungian 

Admi^trative  Contracting  Officer 
Plans  and  Programs  Directorate 
Contracted  Support  Management 


Non-Lincoln  Recipients 

PLEASE  DO  NOT  RETURN 

Permission  is  given  to  destroy  this  document 
when  it  is  no  longer  needed. 


Massachusetts  Institue  of  Technology 
Lincoln  Laboratory 


Preliminary  Design  Review:  GMTI  Narrowband  for  the 
Basic  PCA  Integrated  Radar-Tracker  Application 


A.I.  Reuther 
Group  102 


Project  Report  PCA-IRT-3 

27  February  2003 
Issued  6  February  2004 


Approved  for  public  release;  distribution  is  unUmited . 


Lexington 


Massachusetts 


ABSTRACT 


This  report  describes  ground  moving  target  indicator  (GMTI)  processing  to  be  done  for  the  PCA 
(polymorphism  computing  architecture)  integrated  radar-tracker  application.  GMTI  processing  of  the  raw 
radar  data  is  done  to  extract  targets;  target  reports  are  passed  on  to  the  tracker  part  of  the  integrated  radar- 
tracker  application. 

The  GMTI  processing  described  in  this  report  is  used  to  minimize  the  effects  of  interference  patterns 
and  ground  clutter  and  to  recognize  and  extract  targets  on  the  ground.  The  radar  is  presumed  to  be  using  a 
narrowband  signal.  For  each  data  set,  adaptive  processing  is  used  to  null  interferenee  patterns  and  ground 
clutter  noise.  These  operations  incur  a  high  computational  cost.  Adaptive  processing  requires  the 
calculation  of  an  adaptive  filter  (involving  a  QR  decomposition  and  forward/backsliding)  at  run  time  for 
each  data  cube,  while  target  parameter  estimation  incurs  a  high  computation  cost  as  more  targets  are 
detected  in  each  data  cube. 
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1.  INTRODUCTION 


This  report  describes  ground  moving  target  indicator  (GMTI)  processing  to  be  done  for  the  PCA 
(polymorphism  computing  architecture)  integrated  radar-tracker  application.  GMTI  processing  of  the  raw 
radar  data  is  done  to  extract  targets;  target  reports  are  passed  on  to  the  tracker  part  of  the  integrated  radar- 
tracker  application. 

The  GMTI  processing  described  in  this  report  is  used  to  minimize  the  effects  of  interference  patterns 
and  ground  clutter  and  to  recognize  and  extract  targets  on  the  ground.  The  radar  is  presumed  to  be  using  a 
narrowband  signal.*  For  each  data  set,  adaptive  processing  is  used  to  null  interference  patterns  and  ground 
clutter  noise.  These  operations  incur  a  high  computational  cost.  Adaptive  processing  requires  the 
calculation  of  an  adaptive  filter  (involving  a  QR  decomposition  and  forward/backsliding)  at  run  time  for 
each  data  cube,  while  target  parameter  estimation  incurs  a  high  computation  cost  as  more  targets  are 
detected  in  each  data  cube. 

These  radar  signal  processing  technologies  facilitate  the  extraction  of  targets  from  radar  signal  that 
include  interference  patterns  and  ground  clutter  noise — conditions  in  which  radars  usually  cannot  detect 
any  targets.  Depending  on  the  choices  of  antenna  and  radar  parameters,  most  targets  should  be  detected 
consistently.  However,  some  targets  may  not  be  detected  in  every  detection  interval  because  of  the  trade¬ 
offs  that  exist  in  choosing  these  parameters  and  the  nonlinear  effects  of  adaptive  processing. 

This  report  describes  the  key  features  and  elements  for  a  preliminary  design  review  of  the  GMTI 
narrowband  radar  processing  component  of  the  integrated  moving  target  indicator/tracker  application.  This 
report  includes: 

•  a  high-level  description  of  the  GMTI  narrowband  radar  processing  chain 

•  specifications  of  the  input  and  output  for  each  of  the  GMTI  radar  processing  chain  elements 

•  detailed  descriptions  of  each  of  the  GMTI  radar  processing  chain  elements 

1.1  THE  GMTI  RADAR  PROCESSING  CHAIN 

Before  delving  into  the  details  of  each  of  the  elements  of  the  processing  chain,  a  description  of  the 
entire  GMTI  narrowband  processing  chain  is  in  order.  The  processing  chain  is  illustrated  in  Figure  1.  In 
this  figure,  the  small  numbers  in  the  upper  right  comer  of  the  boxes  refer  to  the  section  in  which  that  stage 
is  discussed  in  greater  detail.  The  shaded  box  with  bold  labels  is  described  in  [1]. 


*  A  wideband  signal  is  a  signal  in  which  the  bandwidth  is  a  significant  fraction  of  the  center  frequency,  while  a 
narrowband  signal  is  one  in  which  the  bandwidth  is  an  insignificant  fraction  of  the  center  frequency.  The  determination 
of  what  constitutes  a  significant  fraction  of  the  center  frequency  is  application  specific. 
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Figure  1.  GMTI  narrowband  radar  processing  chain. 


Before  the  radar  data  is  processed  by  the  GMTI  stages,  a  radar  signal  consisting  of  a  series  of  pulses 
comprising  a  coherent  processing  interval  (CPI)  has  been  transmitted.  The  pulse  repetition  interval  (PRI) 
determines  the  time  interval  between  transmitted  pulses,  and  the  number  of  pulses  per  CPI  is  given  as 
The  signal  reflects  off  of  targets,  the  ground,  and  other  entities  and  is  received  by  the  radar  antenna.  The 
antenna  is  of  subarrays,  each  consisting  of  horizontal  by  Ny  vertical  receiver  elements.  An  array 
of  sensor  inputs,  which  are  also  called  channels,  is  then  derived  from  these  receiver  element  inputs 
using  weighted  summations  of  radiation  energy  measurements.  These  channel  sensor  inputs  are  digitized 
by  A/D  converters  at  a  rate  of  samples  per  PRI.  The  digitization  rate  alludes  to  the  relationships 
between  the  range  dimension  and  the  PRI  dimension:  the  range  dimension  is  referred  to  as  the  fast  time 
dimension,  while  the  PRI  dimension  is  the  slow  time  dimension.  This  x  x  digitized  data  cube 
of  samples  is  the  input  to  the  GMTI  processing  chain. 

The  GMTI  narrowband  processing  chain  consists  of  seven  stages:  time  delay  and  equalization; 
adaptive  beamforming;  pulse  compression;  Doppler  filtering;  space-time  adaptive  processing  (STAP); 
detection;  and  estimation.  This  processing  chain  then  reports  the  resulting  targets  to  the  tracker. 

•  The  time  delay  and  equalization  stage  compensates  for  differences  in  the  transfer  function 
between  channel  sensors. 

•  The  adaptive  beamforming  stage  transforms  the  filtered  data  into  the  beam-space  domain  to 
allow  detection  of  target  signals  coming  from  a  particular  set  of  directions  of  interest  while 
filtering  out  spatially-localized  interference. 
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•  The  pulse  compression  stage  filters  the  data  to  concentrate  the  signal  energy  of  a  relatively  long 
transmitted  radar  pulse  into  a  relatively  short  pulse  response. 

•  The  Doppler  filter  stage  processes  the  data  so  that  the  radial  velocity  of  targets  relative  to  the 
platform  can  be  determined. 

•  The  STAP  stage  is  a  second  beamforming  stage  which  removes  further  interference  and  ground 
clutter  interference. 

•  The  detection  stage  uses  constant  false-alarm  rate  (CFAR)  detection  to  compare  a  radar  signal 
response  to  its  surrounding  signal  responses  to  determine  whether  a  target  is  present  and  uses 
target  grouping  to  eliminate  multiple  target  reports  that  are  actually  just  one  target. 

•  The  estimation  stage  estimates  target  positions  in  order  to  pass  them  to  the  tracking  algorithms. 

•  Finally,  the  target  report  is  passed  to  the  tracking  algorithms  which  are  described  in  [1]. 

1.2  GMTI  CONTROL  INPUT  DATA 

The  GMTI  processing  function  requires  a  control  input  which  is  coded  as  a  Matlab  structure.  The 
names  of  the  elements,  the  variables  they  correspond  to  in  this  report,  and  the  descriptions  of  each  element 
are  given  in  Table  1.  This  structure  must  be  passed  to  the  GMTI  function  when  it  is  called. 
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TABLE  1 

Control  Input  Parameters  for  GMTI  Processing 


Element 

Name 

Description 

Nch 

Number  of  radar  sensor  channels 

Nrg 

Number  of  range  gates  per  radar  data  cube 

Npri 

Npri 

Number  of  pulse  repetition  intervals  per  radar  data 
cube 

^tde 

Ntde 

Number  of  FIR  filter  taps  in  the  time  delay  and 
equalization  filter 

^ab, 

alphaAbf 

Diagonal  loading  factor  in  the  adaptive  beamforming 
stage 

alphaStap 

Diagonal  loading  factor  in  the  STAR  stage 

^pc 

Npc 

Number  of  filter  taps  in  the  pulse  compression  FIR 
filter 

^bm 

Nbm 

Number  of  beams  in  the  radar  data  cube  after  the 
adaptive  beamforming  stage 

Ndop 

Number  of  Doppler  bins  for  Doppler  filtering 

^stag 

Nstag 

Number  of  PRI  staggers  produced  by  Doppler  filtering 
for  use  in  the  STAR  stage 

^abfTS 

NabfTS 

Number  of  training  samples  for  the  adaptive 
beamforming  stage 

^staoTS 

NstapTS 

Number  of  training  samples  for  the  STAR  stage 

^cnb 

Ncnb 

Number  of  clutter  nulled  beams  in  the  radar  data  cube 
after  the  STAR  stage 

Nf  am 

Number  of  receiver  subarrays 

Nx 

Number  of  receiver  elements  in  a  subarray  in  the 
horizontal  dimension 

Ny 

Number  of  receiver  elements  in  a  subarray  in  the 
vertical  dimension 

^cfar 

Ncfar 

Number  of  range  gates  used  in  forming  the  noise 
estimate  from  each  side  of  the  cell  under  test 

G 

G 

Number  of  guard  range  cells  on  both  sides  of  the  cell 
under  test 

mu 

Normalized  target  power  threshold 

The  control  input  parameters  presented  in  the  table  above  are  also  shown  in  tables  at  the  end  of  the 
stage  section  in  which  they  are  used. 
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Beyond  the  control  input  parameters,  there  are  a  number  of  data  product  inputs,  as  shown  in  Table  2. 


TABLE  2 

Data  Product  Inputs  for  GMTI  Processing 


Parameter 

Element 

Name 

Description 

'^bm 

Tbm 

Vector  of  range  gate  indices  for  adaptive  beamforming 
stage  training  samples 

Kc 

hpc 

Coefficient  vector  of  the  pulse  compression  filter 

T 

stap 

Tstap 

Vector  of  space-time  snapshot  indices  for  the  training 
data  for  space-time  adaptive  processing 

1.3  GMTI  INPUT  DATA 

The  input  data  to  the  GMTI  radar  processing  chain  is  a  series  of  three-dimensional  radar  data  cubes 
whose  dimensions  are  channel,  range,  and  pulses  (see  Figure  2).  The  collection  of  the  samples  of  the  radar 
data  cubes  is  explained  in  Section  1.1.  Individual  elements  of  the  radar  data  cube  C  are  referred  to  as  C(i,j, 
k),  where  i  is  the  channel  index, y  is  the  range  index,  and  k  is  the  pulse  index. 

The  series  of  radar  data  cubes  will  be  indexed  using  a  cell  array  in  Matlab:^ 
input_data{ 1 }  =  cubel; 
input_data{2 }  =  cube 2 ; 
input_data { 3 }  =  cube3 ; 
etc . 

Depending  on  the  size  of  the  data  cubes,  they  can  be  generated  separately  and  stored  or  they  can  be 
generated  as  they  are  needed. 


'y 

Items  in  Courier  font  correspond  to  functions  or  variables  in  the  Matlab  code. 
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Njj,  Channels 


Figure  2.  Radar  data  cube,  the  input  data  to  the  GMTl processing. 


1.4  GMTI  OUTPUT  DATA 

The  tracker  receives  the  following  data  for  each  target: 

•  target  power  or  signal-to-noise  ratio  (SNR) 

•  target  azimuth  (in  degrees  or  radians) 

•  target  range  (in  meters  or  kilometers) 

•  target  radial  velocity  (in  m/s  or  km/h) 

•  target  time  stamp 

The  tracker  is  described  in  [1], 
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2.  FUNCTIONAL  DESCRIPTION 


In  this  section,  each  of  the  processing  stages  are  described  in  greater  detail.  For  convenience.  Figure 
3  shows  the  GMTI  narrowband  processing  chain  with  the  sizes  of  the  passed  parameters. 
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Figure  3.  GMTI  narrowband  processing  chain  with  passed  parameter  sizes. 


2.1  TIME  DELAY  AND  EQUALIZATION 

Input:  Data  cubes  of  size  x  x  . 

Output:  Data  cubes  of  size  x  x  . 

The  time  delay  and  equalization  stage  is  applied  to  each  data  band,  and  it  uses  a  finite  impulse 
response  (FIR)  filter  on  each  range  vector  for  each  channel  and  each  pulse  to  compensate  for  signal 
differences  between  channel  sensors.  That  is,  in  Matlab  notation,  it  operates  on  each  x  =  C(i,  k,  s),  where 
i  is  the  channel  index  and  k  is  the  PRI  index.  In  its  simplest  form,  a  single  tap  adjustment  can  be  folded  into 
the  digital  beamforming  weights;  however,  more  complex  equalization  requires  using  programmable  FIR 
filters.  The  FIR  filter  has  filter  coefficients,  ,  0  <  /  <  -  1 .  The  output  of  the  filter,  the  vector  y , 

is  the  convolution  of  the  filter  with  the  input  x : 
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ym  =  E  fori  =  0- -Kg 

1  =  0 


Depending  on  the  number  of  filter  coefficients,  this  filter  can  be  implemented  using  fast 
convolution  with  FFTs  [4]  or  by  explicitly  convolving  the  input  signal  and  the  filter  coefficients  ^ 
depicted  in  the  equation  above. 

The  time  delay  and  equalization  stage  parameter  is  listed  in  Table  3. 


TABLE  3 

Time  Delay  and  Equalization  Stage  Parameter 


Name 

Description 

Kde 

Number  of  taps  in  the  time  delay  and  equalization  fiiters 

2.2  ADAPTIVE  BEAMFORMING 


Input:  Data  cubes  of  size  x  x  . 

Output:  Data  cubes  of  size  x  x  . 

The  adaptive  beamforming  stage  transforms  the  filtered  data  into  the  beam-space  domain  to  allow 
detection  of  target  signals  coming  from  a  particular  set  of  directions  of  interest  while  filtering  out  jamming 
(spatially-localized)  interference.  The  filtering  is  executed  for  each  radar  data  cube  within  each  CPI. 


First,  a  steering  matrix,  FI,  and  antenna  subarray  response  matrix,  Zl,  must  be  determined  for  the 
relative  direction  from  which  the  radar  signal  is  being  received.  For  this  simulation,  the  steering  matrix 
only  involves  the  azimuth  angle  and  ignores  the  elevation  angle  by  assuming  that  it  is  set  to  zero.  The 
vector  of  the  steering  azimuth  angles  is  determined  with  a  linear  distribution  over  the  beam  range  of  the 
radar  with  endpoints  at  the  -3  dB  signal  attenuation  points  of  a  radar  beam.  Assuming  that  the  -3  dB  signal 
attenuation  points  of  a  radar  beam  are  at  ±1°  =  ±7r/180  rad  and  the  look  angle  of  the  radar  is  0,  then  the 


0-  — ...0  +  -^ 

,  .  ,  180  180 
evenly  spaced  angles,  #  has  a  total  of  points. 


vector  of  the  steering  azimuth  angle  is  O  = 


such  that  the  points  in  between  are 


Tl  is  the  antenna  subarray  response  matrix  and  is  of  dimension  where  is  the 

number  of  receiver  elements  (or  receiver  modules)  in  a  subarray  in  horizontal  direction,  and  Nj-^^  are  the 
number  of  horizontal-direction  subarrays.  To  acquire  the  Tl  matrix,  one  must  process  the  antenna 
composition  matrix,  T,  so  that  it  represents  the  received  energy  of  the  look  direction.  The  contents  of  the  T 
matrix  depend  on  the  characteristics  of  the  antenna;  each  channel  is  made  up  of  a  weighted  summation  of 
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the  received  energy  at  each  element.  Hence,  the  elements  in  each  column  of  the  matrix  convey  the 
weighting  of  the  contribution  element  in  the  subarrays  to  that  channel  signal. 


The  T  matrix  must  be  processed  to  reflect  both  the  look  direction  and  the  vector  of  steering  azimuth 
angles.  First,  each  element  of  T  is  processed  for  the  look  direction  such  that  '' 


Tm,n^ 


where  m  =  1 , 


■  ^fam 


)  and  n 


.■Ni 

'  A  '  t,  A,  i  A. 

element  of  h,  and  is  the  Ath  element  of  O. 


.  Then  defining  h  as  a  Hamming  window  vector  of  length 


,  where  /  =  1, ...,  (N^  ■  ,k=l, ...,  Ni,„,  hj  is  the  /th 


FI'  can  now  be  computed  as  FI'  =  (Tl^-ri)  ■  T{\^ ■  U) .  Since  taking  the  inverse  of 
(Tl  ■  T\)  can  be  numerically  unstable,  the  method  described  in  Appendix  A  should  be  employed. 
Defining  =  T\^ ■  {7^,  for  A:  =  1,  ...,  Nf^^,  let  V  =  lqhouse(T\^);  V  is  a  lower  triangular 
matrix.  Now  using  tj^  as  a  temporaiy  vector  of  size  by  1,  we  can  solve  L'tj^  =  for  using  forward- 
substitution  and  (Z')  for  v'^  using  back-substitution.  Then,  FI'  =  [v'j  ...v'^...v'jy  ]  .Finally, 

FI  is  formed  by  taking  the  first  columns  of  the  resultant  Q  matrix  of  the  QR  decomposition  of  FI ' , 
thereby  making  the  columns  of  FI  orthogonal  to  each  other. 


To  perform  adaptive  beamforming,  a  clutter-free  data  matrix  is  extracted  from  C.  Typically, 
the  first  several  range  gates  of  a  CPI  are  considered  clutter-free  so  only  the  interference  signals  are  sensed 
by  the  radar  receiver.  A  total  of  N^tfrs  column  vectors  are  needed,  and  the  choice  of  which  range  samples 
to  use  is  given  by  the  A/'a^y75-length  index  vector  =  {pj, ...,  PMabfrs)-  Commonly,  SN^j^.  Let  x,- 

be  the  rth  column  vector  of  x\  is  formed  by  taking  the  x  matrix  at  range  gate  p,;  in  Matlab 
f^pn  ] 


terms: 


yi  = 


Y  C(:,  Pi,  k) 

U  =  1  y 


which  becomes  a  column  vector  of  length  Nch-  Letting 


6^  =  ^  the  Ath  element  of  the  -length  column  vector  5,  where  A^  =  1, ...,  Np^i,  and  d  is 

^^pri 

the  normalized  Doppler  filter  (a  scalar  number)  that  is  opposite  clutter  which  is  dependent  on  look  angle 
and  system  parameters.  Also,  let  t  be  an  -length  Chebychev  window  filter  column  vector  with  a  desired 

sldelobe  attenuation,  for  instance,  60  dB.  Now  x'^  becomes  the  element-wise  multiplication  ofy,,  6,  and  t 

(Xy-  =  yj.*5.*t,  in  Matlab  terms). 


Now  is  calculated  with  X^^y  such  that 


^abf 


zX 


bfrs 


(2) 
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M  ^ 

and  ^  .  Thus,  an  N {N + N matrix  and  R  is  an  iV^^xiV^^  matrix. 

Note  that  a  diagonal  loading  factor  is  used  to  augment  the  diagonal  of  so  that  it  has  a  lower 

probabiity  of  being  singular. 

To  compute  the  columns  of  the  beamforming,  interference-nulling  matrix,  iri„(ofsizeiV^^xiV^^), 
one  must  solve: 


w 


m 


A— J 

R 


(3) 


where  is  the  mth  column  of  R  is  the  estimated  covariance  matrix,  and  is  the  wth  column  of  the 
steering  matrix,  VI,  computed  above.  The  denominator  of  Equation  (3)  is  a  scalar  normalization  factor. 
The  above  equation  is  referred  to  as  a  Wiener-Hopf  equation  for  calculating  adaptive  filter  weights. 


While  this  explains  the  mathematics  behind  calculating  the  interference-nulling  matrix,  calculating 
the  inverse  of  R  is  not  used  due  to  numerical  stability  issues.  Instead  LQ-decomposition  and  forward-  and 
backward-substitution  is  used  to  calculate  each  of  the  vectors.  As  explained  in  Appendix  A,  an  LQ- 
decomposition  is  performed  on  an  x  matrix:  V  =  lqhouse(A^j^p  in  Matlab; 

V  is  a  lower  triangular  x  matrix.  Now  using  as  a  temporary  vector  of  size  x  1 ,  we  can 

solve  L't^  =  for  using  forward-substitution  and  (L')  w'^  =  for  w'^  using  back-substitution. 
Finally,  compute  account  the  denominator  scalar  normalization  factor. 

Refer  to  Appendix  A  for  a  deeper  explanation  of  this  method  of  solving  this  Wiener-Hopf  equation. 

Finally,  is  used  to  beamform  the  radar  data  cube.  Using  Matlab  notation  in  defining  the  matrix 
X^,  =  C(:, k)  with  size  x  ,  the  matrix  7^  as  the  output  of  the  beamforming, 


(4) 


Thus,  the  matrix  7^  is  of  size  ^  •’  “  ^k-  Figure  4  shows  a  conceputal  data 

cube  and  its  partitioning  for  this  operation. 
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Figure  4.  Data  partitioning  of  input  data  cube  for  adaptive  beamforming  stage. 


The  data  parameters  of  the  adaptive  beamforming  stage  are  listed  in  Table  4. 

TABLE  4 

Adaptive  Beamforming  Stage  Parameters 


Name 

Description 

^abf 

Adaptive  beamforming  diagonal  loading  factor 

^abfrS 

Number  of  training  data  set  vectors  in  adaptive  beamforming  stage 

'^bm 

Vector  of  range  gate  indices  for  adaptive  beamforming  stage  training  samples 

2.3  PULSE  COMPRESSION 

Input:  Data  cubes  of  size  x  x  . 

Output:  Data  cubes  of  size  x  _  j )  x  . 

The  pulse  compression  stage  filters  the  data  to  concentrate  the  signal  energy  of  a  relatively  long 
transmitted  radar  pulse  into  a  relatively  short  pulse  response.  This  compression  is  accomplished  by 
applying  a  FIR  filter  operation  to  the  range  data  for  each  pulse  and  channel  within  each  radar  data  cube. 
That  is,  in  Matlab  notation,  it  operates  on  each  x  =  C(i,  k),  each  of  which  are  of  length  this 
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partitioning  is  illustrated  in  Figure  5,  The  FIR  filter  has  coefficients  /e  1}.  The 

output  of  the  filter,  the  vector  y  of  length  +  Np^  -  1,  is  the  convolution  of  the  filter  with  the  input 
x: 


T[/]  =  X  x\j-l]hp^[l],forj  = 
1  =  0 


(5) 


Figure  5.  Data  partitioning  of  input  data  cube  for  pulse  compression  stage. 


Unlike  the  time  delay  and  equalization  stage,  the  number  of  taps  in  the  pulse  compression  FIR  filter 
is  usually  much  larger  than  four  or  so  taps.  Therefore,  it  is  usually  more  efficient  to  implement  this 
convolution  with  the  fast  convolution  technique  using  FFTs,  possibly  even  considering  an  overlap-and- 
save  method  [4], 

The  pulse  compression  stage  data  parameter  is  listed  in  Table  5. 

TABLE  5 


Pulse  Compression  Stage  Parameter 


Name 

Description 

^pc 

Number  of  taps  in  the  subband  pulse  compression  filter 
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2.4  DOPPLER  FILTERING 


Input:  Data  cubes  of  size  x  -  1 )  x  . 

Output:  Data  cubes  of  size  x  -  1 )  x  x  .  In  the  Matlab  implementation, 

the  subband  data  cube  is  coded  to  be  of  size  x  -  1 )  x  . 

The  Doppler  filter  stage  processes  the  data  so  that  the  radial  velocity  of  targets  relative  to  the 
platform  can  be  determined.  Doppler  filtering  accomplishes  this  using  a  straightforward  FFT  applied  to  the 
pulse  data  for  each  range  gate  and  beam.  A  further  step  made  in  the  Doppler  filtering  which  aids  the  STAP 
stage  in  filtering  out  ground  clutter  involves  taking  staggered  pulse  sample  sets;  is  the  number  of 
pulse  staggers.  Using  g  (1  <g<  )  as  the  index  of  the  stagger  number,  the  staggers  generated  in 

Matlab  notation  are  y^g  =  C(i,  j,  g-'^pri-  ^stag  +  g))'->  yij,g  a  vector  of  length  .  For  example,  with 
^stag  ^  ^  input  stagger  vector  would  include  the  first  pulse  sample  up  to  the  Ap^^.-2nd  pulse 

sample,  the  second  input  stagger  vector  would  include  the  second  pulse  sample  up  to  the  -  1  st  pulse 
sample,  and  the  third  input  stagger  vector  would  include  the  third  pulse  sample  up  to  the  A^^^  th  (last) 
pulse  sample.  This  is  illustrated  in  Figure  6.  With  A^^^g  =  1,  there  is  only  one  stagger  which  is  the  entire 
set  of  pulse  samples  for  a  given  beam  and  range  gate. 

To  smooth  out  the  ringing  effect  of  the  FFT,  a  taper  is  applied  to  each  g  vector.  Let  t  be  an  Np^,- 
length  Chebychev  window  filter  column  vector  with  a  desired  sidelobe  attenuation,  for  instance,  60  dB. 
The  tapered  y,j  g  vector  is  an  element-wise  product:  y'- j  g-  (j,-  j  g)  The  output  of  the  Doppler  filter 
is 


''.A? 


FFTiy'.jJ 


(6) 


for  each  beam,  range  gate,  and  stagger.  Each  y  g  is  a  vector  of  size  A^^^ .  These  vectors  can  be  stored 
in  a  variety  of  ways;  in  the  Matlab  implementation,  the  Doppler  filtered  staggers  are  indexed  using  the 
channel/beam  dimension.  That  is,  the  index  of  the  first  dimension  of  C{i,  j,  k)  accesses  all  of  the  beams  of 
the  first  stagger  (g  =  1,  /  =  1, ...,  followed  by  all  of  the  beams  of  the  second  stagger  (g  =2,  i  =  + 

1, ...,  2A^^)  and  so  on. 

The  parameters  of  the  Doppler  processing  stage  are  listed  in  Table  6. 
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ll'^^er  ^  -  samrtey  2  »ru 

1 

1  ^ta^er  3'thnj 

Figure  6.  Data  partitioning  shewing  three  staggers  of  input  data  cube  for  Doppler  filtering  stage. 


TABLE  6 

Doppler  Processing  Stage  Parameters 


Name 

Description 

M 

^^stag 

Number  of  pulse  staggers 

2.5  SPACE-TIME  ADAPTIVE  PROCESSING  (STAP) 


Input:  Data  cubes  of  size  x  -  1 )  X  x  .  In  the  Matlab  implementation, 

the  subband  data  cube  is  coded  to  be  of  size  ( x  -  1 )  x  . 


Output:  Data  cubes  of  size  x  -  1 )  x  . 


The  STAP  stage  is  a  second  beamforming  stage  which  removes  further  radar  spatially-localized 
interference  and  ground  clutter  (space-  and  time-distributed)  interference.  Using  the  z.  j  ^  outputs  from 
the  Doppler  filtering  stage,  define 


^J,k,g 


(7) 


14 


In  other  words,  Xjj^g  is  constituted  with  the  Nf,^  complex  elements  for  a  given  stagger,  range  gate, 
and  Doppler  gate.  The  size  of  Xj  /t  g  ^bm  ^  ^  stagger  of  samples.  Next,  the  staggers  are 

stacked  upon  each  other  to  get  a  space-time  snapshot: 


j,  k,  +  Ij 


(8) 


for  the  test  cells,  which  makes  Xj  ^  an  ^bm^stag  ^  ^  vector.  We  will  use  these  vectors  of  size 
^bm^stag  ^  1  our  input  data  to  the  STAP  fdter,  but  we  shall  first  build  STAP  steering  vectors,  and  then 
compute  the  adaptive  filter  matrix. 

As  in  the  adaptive  beamforming  stage,  a  steering  matrix,  VI,  and  antenna  subarray  response  matrix, 
Tl,  must  be  determined  for  the  relative  direction  from  which  the  radar  signal  has  been  received.  J2  is 
calculated  from  the  matrix  multiplication  of  T\  and  Wj^^\  T2  =  T\  ■  ■  T2  is  of  dimensions 

'  ^fan^  ^  ^bm  ’  "^^ere  Ny.  is  the  number  of  receiver  elements  in  a  subarray  in  horizontal  direction,  and 
are  the  number  of  horizontal-direction  subarrays. 

Also  as  in  the  adaptive  beamforming  stage,  the  vector  of  steering  azimuth  angles  is  determined  with 
a  linear  distribution  over  the  beam  range  of  the  radar  with  endpoints  at  the  -3  dB  signal  attenuation  points 
of  a  radar  beam.  Assuming  that  the  -3  dB  signal  attenuation  points  of  a  radar  beam  are  at  ±1  °  =  ±7i/l  80  rad, 
and  the  look  angle  of  the  radar  is  0,  then  the  vector  of  steering  azimuth  angles  is 
^  ~  ^  ~  Tio  ^  points  in  between  are  evenly  spaced  angles.  But  unlike  the 

adaptiW  beamforming  stage,  O  has  a  total  of  points,  instead  of  points.  Then,  as  in  the  adaptive 
beamforming  stage,  defining  /i  as  a  Hamming  window  vector  of  length  {N  ■  Nr  ^ ,  we  can  define 

““/Th/c))  ^ 

Ui  ^  =  h<S>i^e  ,  where  /  =  1, ...,  {N^  •  ,  A  =  1, ...,  /i,  is  the  zth  element  of  h,  and  is  the 

Ath  element  of  O. 

V2'  can  now  be  computed  as 


V2'  =  {T2^  ■  T2)  *  •  {T2^  ■  U)  .  (9) 

Since  taking  the  inverse  of  (r2  -72)  can  be  numerically  unstable,  the  method  described  in 
Appendix  A  should  be  employed.  Defining  V]^  =  T2  ■  Uj^,  for  k  =  1,  ...,  N^„b’ 

L'=  lqhouse{T2  );  L'  is  a  lower  triangular  matrix.  Now  using  as  a  temporary  vector  of 

size  X  1 ,  we  can  solve  L'tj^  =  for  using  forward-substitution  and  {L')  v'^  =  tj^  for  v'^  using 
back-substitution.  Then,  F2'  =  [v'j ...v'^...v'^  ]. 
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To  this  point,  all  of  the  processing  applies  to  the  subband  data  cube  as  a  whole,  but  the  rest  of  the 
processing  must  be  completed  for  each  Doppler  bln,  8.  Within  each  Doppler  bin  to  accommodate  the  STAP 
calculation  for  each  stagger,  V2'  is  replicated  with  an  added  phasor  term  which  is  dependent  on  the 
Doppler  bin,  8.  Letting  the  phasor  equal: 


r 


(|)(g)  =  exp 


g- 


N. 


dop 


modi 

JJ 


J 


(10) 


where  g=\, ...,  V2"  is  calculated  as: 


F2"  =  [^(1)  •  F2'...^(g)  •  V2'..MNgtag^  ■  F2']  . 


(11) 


F2"  is  a  x  matrix.  Finally,  V2  is  formed  by  taking  the  first  columns  of  the 

resultant  Q  matrix  of  the  QR  decomposition  of  F2" ,  thereby  making  the  columns  of  F2  orthogonal  to  each 
other.  F2  is  a  x  matrix.  The  first  rows  of  the  R  matrix  of  the  QR  decomposition  of 

F2'  are  also  retained  and  renamed  R2,  which  is  a  x  matrix. 

For  each  Doppler  bin,  a  sample  data  matrix  is  extracted  from  the  output  of  the  Doppler 

filtering  stage.  A  total  of  Ngi^pps  column  vectors  are  needed,  and  the  choice  of  which  range  samples  to  use 
is  given  by  the  AT^^pT-s-length  index  vector  =  {pj, ...,  } .  Commonly,  Ngtapjs  =  ^^bnf^stag- 

Using  the  xy  ^  vectors  of  size  x  1 ,  which  presented  at  tiie%eginning  of  this  section 


X. 


(8)  _ 


stop 


- 


Pat,  S. 

stap 


(12) 


Now  Agfg^^^  is  calculated  with  XgfgJ-^\  such  that 


stap 


(6) 


stop 


J^stapTS 


X. 


(5) 


stapTS 


J' 


^stap^ 


(13) 


'(8) 


and  R^'"  =  AgJ^^Ag^g^^Y-  Thus,  is  an  N g^^^^  {N N g^^^)  matrix,  and  M''' 

is  an  N/^g^Ng^gg  x  matrix.  Note  that  a  diagonal  loading  factor  Ugigj^  is  used  to  augment  the 

diagonal  of  so  that  it  has  a  lower  probability  of  being  singular. 

To  compute  the  columns  of  the  beamforming,  interference-nulling  matrix,  (of  size 

^bm^stag  ^  ^cnb)^  o"®  »"ust  solve: 


(8) 


H 


W 


8,6 


''(5)  “1 

(R  )  ^6 

Vfc  (R  ) 


16 


(8) 

where  wg  ^  is  the  feth  column  of  R  is  the  estimated  covariance  matrix,  and  is  the  bi\i  column 
{b=\,  of  the  steering  matrix,  V2,  computed  above.  The  denominator  of  Equation  (14)  is  a  scalar 

normalization  factor.  The  above  equation  is  referred  to  as  a  Wiener-Hopf  equation  for  calculating  adaptive 
filter  weights. 


While ^^is  explains  the  mathematics  behind  calculating  the  clutter-nulling  matrix,  calculating  the 
inverse  of  R  is  not  used  due  to  numerical  stability  issues.  Instead,  LQ-decomposition  and  forward-  and 
back-substitution  is  used  to  calculate  each  of  the  w^j,  vectors  (b  =  \,  ...,  N^„i,)  of  the  matrix.  As 

explained  in  Appendix  A,  an  LQ-decomposition  is  performed  on  ^stap^’ 


^ brn^ stag  X  ^bm 

H,bh  " 


stapTS  ^  ^bm^stag^ 


matrix:  Egg  =  lqhouse(^^,ap*^^^;  T, 


stop 


matrix.  Now  using  t/,  as  a  temporary  vector  of  size  x  1 

Slug  ^  H  ^ 


:.g  g  is  a  lower  triangular 
we  can  solve 


for  tjj  using  forward-substitution  and  (Zg  g)  w'g  ^  = 


bm  stag ' 

for  w'g  ^  using  back-substitution. 


factor. 


Lastly,  compute  Wg  ^  =  w'g  ^/((tg)  to  take  into  account  the  denominator  scalar  normalization 


Refer  to  Appendix  A  for  a  more  detailed  explanation  of  this  method  of  solving  the  Wiener-Hopf 
equation. 

Finally,  the  adaptive  weights  vectors  are  applied  to  the  snapshots  Xj  g  of  size  X  1- 


5 


(15) 


6  =  I...N 


The  resulting  vector  yj  g  is  of  size  x  1  and  must  be  calculated  for  each  j  =  1  and 


dop  ' 


The  STAP  stage  data  parameters  are  listed  in  Table  7. 


TABLE  7 

Space-Time  Adaptive  Processing  Stage  Parameters 


Name 

Description 

M 

stag 

Number  of  pulse  staggers 

^stapTS 

Number  of  training  data  set  vectors  in  the  STAP  stage 

'^stap 

Vector  of  space-time  snapshot  indices  for  the  training  data  set  for  space-time 
adaptive  processing 
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2.6  TARGET  DETECTION 


Input:  Processed  radar  data  cube  of  size  x  (N^g  +  -  1 )  x  . 

Output:  Coordinates  of  detected  targets  in  radar  data  cube. 

Though  there  are  many  techniques  for  detecting  targets  in  the  processed  radar  data  cube  including 
various  integrators  and  hard  limiters,  this  application  uses  a  version  of  constant  false  alarm  rate  (CFAR) 
detection  with  three-dimensional  grouping.  The  data  parameters  of  the  target  detection  stage  are  listed  in 
Table  8. 


TABLE  8 


Target  Detection  Stage  Parameters 


Name 

Description 

^cfar 

Number  of  range  gates  used  in  forming  the  noise  estimate  from  each  side  of 
the  cell  under  test 

G 

Number  of  guard  cells  on  both  sides  of  the  cell  under  test 

P 

Normalized  target  power  threshold 

2.6.1  CFAR  Detection 

During  CFAR  detection,  a  local  noise  estimate  is  computed  from  the  range  gates  near  the  cell 

C(i,  j,  k)  under  test.  A  number  of  guard  gates  G  immediately  next  to  the  cell  under  test  will  not  be  included 
in  the  local  noise  estimate.  (This  number  does  not  affect  the  throughput.)  The  range  cells  involved  in 
calculating  the  noise  estimate  for  a  particular  Doppler  bin  and  beam  are  shown  in  Figure  7.  For  each  cell 
C(i,  j,  k),  the  value  of  the  noise  estimate  T(i,  j,  k)  is  calculated  as: 

T(i,j,k)  =  ,^  X  \C{i,j  +  l,kf  +  \C{i,j-l,kf  .  (16) 

At  the  boundaries  when  j-{G  +  -  1 )  <  1  or  j  +  (G  +  >  N,  only  the  side  range  gate 

cells  that  are  within  the  data  cube  are  included,  and  the  denominator  of  the  leading  scaling  factor  is 
decreased  to  reflect  the  number  of  side  range  gate  cells  that  were  actually  used  for  the  calculation.  For  each 
cell  C(i,  J,  k),  the  quantity  \C(i,j,  k)^/T{i,j,  k)  is  calculated;  this  represents  the  normalized  power  in  the 
cell  under  test.  If  this  normalized  power  exceeds  a  threshold  p ,  the  cell  is  considered  to  contain  a  target. 
This  threshold  can  be  adaptively  determined,  but  for  this  application,  it  is  set  to  a  constant  The  only  cells 
in  the  data  cube  that  are  processed  by  the  following  stages  are  those  where  a  target  h^  been  detected. 
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Figure  7.  Data  partitioning  of  input  for  CFAR  detection  stage. 


2.6.2  Three-dimensional  Grouping 

During  this  step,  detection  in  adjacent  cells  are  grouped  to  avoid  having  multiple  detections 
associated  with  the  same  target.  The  power  of  each  target  detected  by  the  CFAR  algorithm  will  be 
compared  to  the  power  of  each  cell  in  the  3x3x3  cube  centered  on  that  cell.  Therefore,  each  cell  under 
test  will  be  compared  to  26  other  cells  (27  -  1 ,  where  27  =  3  and  1  cell  is  the  cell  under  test,  the  cell  in 
the  middle  of  the  cube).  If  the  detection’s  raw  power  is  the  local  maximum,  it  is  retained.  Otherwise,  the 
detection  is  discarded;  it  is  considered  grouped  with  the  detection  in  the  cell  with  the  higher  power. 

2.7  TARGET  PARAMETER  ESTIMATION 

Input:  Coordinates  of  detected  targets  in  radar  data  cube. 

Output:  A  target  report  for  each  detected  target  with  sub-bin  resolution. 

This  stage  of  the  GMTI  processing  packages  the  target  parameters  for  the  target  report.  As  presented 
in  Section  1 .4,  a  target  report  consists  of  the  following  data  for  each  target: 

•  target  power  or  signal-to-noise  ratio  (SNR) 

•  target  azimuth  (in  degrees  or  radians) 

•  target  range  (in  meters  or  kilometers) 
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target  radial  velocity  (in  m/s  or  km/h) 
target  time  stamp 


2 

The  target  power  is  the  quantity  \C{i,j,k)\  /T(i,J,k)  for  each  target  cell  C(i,  j,  k)  which  was 
calculated  by  the  CFAR  detection.  The  target  time  stamp  will  be  reported  as  the  input  time  stamp  of  the 
data  cube  from  which  this  target  was  determined. 

For  the  target  azimuth,  range,  and  radial  velocity,  sub-bin  estimates  can  be  determined  by 
interpolating  between  the  processed  signal  return  points  near  the  target.  These  interpolation  estimators 
require  a  vector  of  data  between  which  these  interpolations  will  be  calculated.  A  vector  of  interpolation 
points  over  the  vector  span  are  determined,  and  from  these  interpolation  points,  the  sub-bin  target  location 
is  where  the  maximum  value  is  found  in  this  interpolation  vector.  For  this  simulation  software,  the  target 
azimuth  angle  estimate  is  conducted  with  a  maximum  likelihood  estimator,  while  the  target  range  and 
radial  velocity  (Doppler)  estimates  are  conducted  with  a  cubic  spline  interpolation  estimator.  The 
following  two  sections  describe  these  two  estimators.  Finally,  these  sub-bin  estimates  are  converted  into 
their  corresponding  target  data  values. 

2.7.1  Maximum  Likelihood  Estimator 

To  perform  maximum  likelihood  estimation  (MLE)  of  the  azimuth  angle,  post-Doppler  filtering,  pre- 
STAP  data  must  be  used.  Two  hundred  evenly-spaced  azimuth  angle  samples  are  taken  across  the  entire 
azimuth  range,  and  the  index  of  the  angle  at  which  the  maximum  value  occurs  is  returned  as  the  sub-bin 
estimate.  All  of  these  calculations  must  be  made  for  each  target  report  generated  by  the  detection  range. 

The  MLE  function  takes  as  input  a  vector  of  azimuth  values,  xjj^  an  x  1  vector  as 

calculated  in  Equation  (8),  where  the  target  in  question  has  been  detected  in  range  bin  j  and  Doppler  bin  k. 
It  also  takes  a  vector  of  indices  (azimuth  bin  numbers  which  correspond  to  the  vector  of 

azimuth  angles.  (They  must  be  of  the  same  length.) 

Similar  to  the  STAP  stage,  a  steering  matrix,  V2,  and  antenna  subarray  response  matrix,  72,  must  be 
determined  for  the  relative  direction  from  which  the  radar  signal  is  being  received.  72  is  calculated  from 
the  matrix  multiplication  of  T\  and  T2  =  Tl  ■  72  is  of  dimensions  (A^  •  A^^^)  x  , 

where  is  the  number  of  receiver  elements  in  a  subarray  in  horizontal  direction,  and  are  the  number 
of  horizontal-direction  subarrays. 

Also  similar  to  the  STAP  stage,  the  vector  of  steering  azimuth  angles  is  determined  with  a  linear 
distribution  over  the  beam  range  of  the  radar  with  endpoints  at  the  -3  dB  signal  attenuation  points  of  a 
radar  beam.  Assuming  that  the  -3  dB  signal  attenuation  points  of  a  radar  beam  are  at  ±1  °  =  ±Jt/l  80  rad,  and 
the  look  angle  of  the  radar  is  0,  then  the  vector  of  steering  azimuth  angles  is  . . .  0  +  ■—  , 

such  that  the  points  in  between  are  evenly  spaced  angles.  But  unlike  the  STAP  stage,  ®  has  a  total  of  2u0 
points,  instead  of  N^„i,  points.  Then,  as  in  the  STAP  stage,  defining  h  as  a  Hamming  window  vector  of 
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length  {N^  ■  ,  we  can  define  t/.^ ^  ,  where  i  =  1, ...,  {N^  ■  ,k=\, ...,  200,  /i,-  is 

the  /th  element  of  h,  and  is  the  ^h  element  of  O. 

V2'  can  now  be  computed  as 


V2'  =  {T2^ ■  T2)  *  •  {T2^ ■  U) 


(17) 


Since  taking  the  inverse  of  (T2  ■  T2)  can  be  numerically  unstable,  the  method  described  in 
Appendix  A  should  be  employed.  Defining  vj^  =  T2  ■  Uj^,  for  k  =  \,  ..., 

L'=  lqhouse{T2  );  L'  is  a  lower  triangular  200  matrix.  Now  using  as  a  temporary  vector  of 

size  A^^  X  1 ,  we  can  solve  L'tj^  =  Vj^  for  using  forward-substitution  and  (/,')  v'^  =  tj^  for  v'^  using 
back-substitution.  Then,  V2'  =  [v',  ...v'.  ...v'^  ]. 

Since  the  target  detection  occurred  in  a  specific  Doppler  bin,  these  calculations  must  only  be 
conducted  for  that  specific  Doppler  bin,  8.  Within  the  calculations  for  that  Doppler  bin,  V2'  is  replicated 
with  an  added  phasor  term  which  is  dependent  on  the  Doppler  bin,  8.  Letting  the  phasor  equal: 


(t)(g)  =  exp 


g- 


vv 


A 


dop 


modi 


JJ 


(18) 


where  g-=  1, ...,  V2  is  calculated  as: 


F2  =  [^(D-  F2'...4)(g).  F2'...(t)(A^,^p-  F2'] 


(19) 


F2  is  a  (A^^  •  A^^^^)  x  200  matrix.  Using  the  A^f^p  matrix  for  Doppler  k  calculated  in  Equation 

(13),  we  can  find  Wq  by  solving  the  Wiener-Hopf  equation:  A^i^p^^^Wq  =  F2  for  Wq.  Refer  to  Appendix  A 
and  paragraphs  above  to  formulate  this  in  a  computationally  stable  manner. 

Wq  is  then  used  in  the  MLE  equation  to  find  a  vector  of  sub-bin  azimuth  angle  return  estimates,  0,  a 
200  element  vector: 


0 


(WQfv2 


(20) 


The  maximum  value  location  across  the  200  elements  of  0  must  then  be  found  and  the 
corresponding  sub-bin  index  value  is  returned. 
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2.7.2  Cubic  Spline  Interpolation  Estimator 


The  cubic  spline  interpolation  estimator,  used  to  determine  sub-bin  range  and  Doppler  estimates, 
takes  as  input  a  vector  of  values  over  which  to  estimate  and  a  vector  indices  (bin  numbers)  which 
correspond  to  the  vector  of  values.  (They  must  be  of  the  same  length.)  These  vectors  are  assembled  by 
taking  each  target  bin  position  and  including  several  points  to  each  side  of  the  target  bin  position  within  the 
dimension  to  be  determined.  So  if  a  sub-bin  range  estimate  is  desired,  and  if  the  detection  stage  has  found 
a  target  in  azimuth  bin  i,  range  bin  p,  and  Doppler  bin  k,  and  we  are  including  three  range  points  to  either 
side  of  the  target,  the  vector  of  values  would  be  jc  =  C(i,  p  -  3:p  +  3,  k),  in  Matlab  terms,  and  the  index 
vector  would  contain  [p  -  3:p  +  3]. 

Cubic  spline  inteipolation  involves  interpolating  with  a  cubic  equation  between  each  set  of  adjacent 
points.  The  technique  constrains  the  interpolation  intervals  to  have  continuous  second  derivatives  which 
makes  the  complete  interpolation  curve  over  all  of  the  value  points  smooth.  The  second  derivatives  can  be 
computed  for  each  of  the  points  in  the  values  vector,  and  the  resulting  system  of  equations  which  solve  for 
the  coefficients  of  other  interpolation  intervals  can  be  assembled  as  a  tridiagonal  matrix  solution.  For  a 
complete  explanation  of  cubic  spline  interpolation,  see  Section  3.3  of  [5]. 

To  evaluation  sub-bin  estimates,  one  must  construct  a  vector  of  sub-bin  indices;  for  instance,  we  can 
subdivide  the  interval  of  indices  into  a  total  of  100  evenly-spaced  points  for  which  we  want  interpolated 
values.  With  the  cubic  interpolation  equation  coefficients  of  each  interval,  we  can  then  evaluate  the 
interpolated  values  at  each  of  the  100  points,  find  the  maximum  value  location,  and  return  the 
corresponding  sub-bin  index  value.  All  of  these  calculations  must  be  made  for  each  target  report  generated 
by  the  detection  stage. 
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3.  MATLAB  IMPLEMENTATION 


The  GMTI  processing  chain  will  be  initialized  with  a  function  with  the  following  signature: 
Gmti_param  =  GMTI_Init(); 

This  function  will  load  control  parameters  outlined  in  Section  1 .2  into  the  structure  GMTI_params. 
The  GMTI  processing  chain  will  be  invoked  with  a  function  with  the  following  signature: 

targets  =  GMTI_run (input_data,  GMTI_params) ; 

This  function  call  passes  the  aforementioned  GMTI_params  and  the  input_data  outlined  in 
Section  1.3,  and  it  returns  the  target  structure  outlined  in  Section  1.4. 

In  order  to  keep  memory  usage  to  a  reasonable  level,  the  GMTI  function  will  process  only  one  radar 
data  cube  at  a  time,  and  within  that  radar  data  cube,  it  will  process  only  one  subband  radar  data  cube  at  a 
time.  Each  of  the  stages  depicted  in  Figure  1  will  have  its  own  function  call  within  the  GMTI  function. 
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APPENDEX  A 

SOLVING  THE  WIENER-HOPE  EQUATION 


Wiener-Hopf  equation  is  used  in  both  the  adaptive  beamforming  and  STAP  stages.  The  equation 
w  =  J?  V  defines  an  adaptive  filter  vector,  w,  as  the  product  of  the  inverted  estimated  noise  covariance 
matrix,  R ,  and  the  nonadaptive  steering  vector,  v.  R  is  composed  using  a  set  of  training  vectors,  Xj.,  such 
thatJ^  =  . 

However,  R  and  its  inverse  are  seldomly  used  to  solve  the  Wiener-Hopf  equation  [3]  because 
calculating  R  and  its  inverse  causes  numerical  instability.  Instead,  LQ-decomposition  with  forward-  and 
back-substitution  is  used.  First,  the  more  general  solution  of  the  STAP  Wiener-Hopf  equation  is  discussed, 
and  then  the  adaptive  beamforming  solution  (in  which  a  diagonal  loading  factor  is  added)  is  shown. 

Compute  the  LQ  decomposition  Xj,  =  LQ,  where  L  becomes  a  lower  triangular  matrix  and  Q 
becomes  a  unitary  matrix.  The  LQ  decomposition  is  related  to  the  QR  decomposition  described  in  [2],  by 
inputting  the  Hermitian  of  Xj,  and  receiving  the  output  as  the  Hermitian  of  L.  Matlab  code  for  the  LQ 
decomposition  follows: 

%lqHouse  Perform  LQ  decomposition  of  input  matrix 
function  L  =  IqHouse (A) 

%  Find  the  size  of  the  input  matrix  A. 

[m,  n]  =  size  (A)  ; 

%  Loop  over  each  row  of  A. 
for  row  =  l:m 


%  Compute  the  Householder  vector  v. 
x=0; 

x(row:n,  1)  =  (A (row,  row:n) ) ' ; 
v  =  x; 

v(row)  =  v(row)  +  (x(row)  /  abs(x(row)))  *  norm(x); 
beta  =  -  2  /  (v'  *  v) ; 


%  Apply  the 
matrix . 

w ( row : m , 
A ( row : m , 


Householder  reflection  to  the  remainder  of  the 

1)  =  beta  *  A{row:m,  row:n)  *  v(row:n); 

row:n)  =  A(row:m,  row:n)  +  (w(row:m)  *  v(row:n)'); 
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end 


for  row 


L  =  A; 

Now  the  estimated  noise  covariance  matrix  can  be  expressed  as: 

R  =  =  (LQ)(LQf  =  LQ(Q^L^)  =  LL^  ,  (21) 


and  the  Wiener-Hopf  equation  can  be  rewritten  as: 


— 1  //  -1 

w  =  R  V  =  {LL  )  V 


{LL^)w  =  (Ll\lL^)  ^  =  V 


(22) 

(23) 


H 

Letting  t  =  L  w,  Lt  =  v  can  be  solved  for  t  by  using  forward-substitution  since  L  is  lower 
triangular.  Then  usi^  t  =  L  w ,  a  solution  for  the  adaptive  filter  vector,  w ,  can  be  found  by  using  back- 
substitution  since  L  is  upper  triangular. 


Now  the  added  complication  of  the  added  diagonal  loading  factor  in  the  adaptive  beamforming  stage 


is  discussed.  To  solve  thi^jproblem,  define  Xj! 


equation,  w  =  ((Z.’)(Z,')  )  v,  can  be  solved  as  described  above. 


ija 


Now  Xj!  =  L'Q ,  and  the  Wiener-Hopf 


Finally,  the  denominator  of  equation  (3)  is  a  scalar  which  can  be  solved  in  the  following  manner. 
Using  R  = 


H — 1  H  H  H  //  -1  -1 

V  R  V  =  v"((L')(^')  )  V  =  V  ((f)  )  (f)  V  . 


|..etting  u  =  (V)  v,  L'u  =  v,  which  can  be  solved  for  u  using  forward-substitution,  and 


^  V  =  u ,  which  is  a  scalar. 
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ACRONYMS 


CFAR  -  Constant  False  Alarm  Rate 
CPI  -  Coherent  Processing  Interval 
FFT  -  Fast  Fourier  Transform 
FIR  -  Finite  Impulse  Response 
GMTI  -  Ground  Moving  Target  Indicator 
PCA  -  Polymorphous  Computer  Architecture 
PRI  -  Pulse  Repetition  Interval 
SNR  -  Signal-to-Noise  Ratio 


STAP  -  Space-Time  Adaptive  Processing 
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